Effect of spin-orbit coupling on the excitation spectrum of Andreev billiards 
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We consider the effect of spin-orbit coupling on the low energy excitation spectrum of an Andreev 
billiard (a quantum dot weakly coupled to a superconductor), using a dynamical numerical model 
(the spin Andreev map). Three effects of spin-orbit coupling are obtained in our simulations: In 
zero magnetic field: (1) the narrowing of the distribution of the excitation gap; (2) the appearance of 
oscillations in the average density of states. In strong magnetic field: (3) the appearance of a peak in 
the average density of states at zero energy. All three effects have been predicted by random-matrix 
theory. 
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I. INTRODUCTION 

A quantum dot in a two-dimensional electron gas has 
a mean level spacing which is independent of energy and 
depends only on geometrical factors (area) and material 
properties (effective mass). The nature of the electron 
dynamics (chaotic versus integrable) and the presence or 
absence of symmetries, such as time-reversal and spin- 
rotation symmetry, have no effect on the mean density of 
states. The situation changes if the quantum dot is cou- 
pled to a superconductor (see Fig. Q}. The presence of 
the superconductor strongly affects the excitation spec- 
trum of such an Andreev billiard. The density of states at 
the Fermi level is suppressed in a way which is sensitive 
to the nature of the dynamics and existing symmetries. 1 
While the effect of broken time-reversal symmetry on the 
density of states has been studied extensively^^ the 
effect of broken spin-rotation symmetry due to spin-orbit 
coupling has only been partially investigated^^ 

In this paper we study in computer simulations three 
effects of the spin-orbit coupling. All of these effects have 
been predicted by random-matrix theory (RMT)j£i2»i£ 
but have so far not been confirmed in a dynamical model. 
The first two of these effects present in the absence 
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Figure 1: Sketch of an Andreev Billiard: A quantum dot (N) 
connected to a superconductor (S) by a point contact. Spin- 
orbit coupling is present in the quantum dot. 



of a magnetic field, are the reduction of the sample-to- 
sample fluctuations of the excitation gap and the appear- 
ance of oscillations as a function of energy in the average 
density of states. The third effect appears in a magnetic 
field strong enough to close the excitation gap. While 
in the absence of spin-orbit coupling the average density 
of states vanishes at the Fermi level, in the presence of 
spin-orbit coupling it peaks at the Fermi level at twice 
the value in the normal stated 



II. PREDICTIONS OF RANDOM-MATRIX 
THEORY 

We begin by briefly summarizing the RMT of the An- 
dreev billiard. 1 In perturbation theory the density of 
states of non-degenerate levels (in zero or weak magnetic 
field) has a square root dependence on energy near the 
gap^iL 
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Here 5 = 2-Kh 2 /mA is the mean level spacing in the 
isolated quantum dot (area A, effective mass m), N = 
Int[fci?IF/7r] is the number of modes in the ballistic point 
contact (width W, Fermi wave vector hp) connecting it 
to the superconductor, and Et — NS/^tt is the Thouless 
energy. These parameters refer to two-fold degenerate 
levels and modes, corresponding to s = 2. If both spin- 
rotation and time-reversal symmetries are broken, the 
two- fold degeneracy is lifted and one should take s = 1. 
The numerical coefficients c and d are magnetic field de- 
pendent. For B = one has 



c = 2 7 5 / 2 , d = (5 - 2y/5) 
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with 7 = (y/5 — l)/2 the reciprocal of the golden ratio. 
For B =/= they should be calculated from the RMT 
solution given in Ref. 0. 

The perturbation theory has S/Et as a small param- 
eter and gives the density of states with an energy reso- 
lution of order Et, which is a macroscopic energy scale. 
Since spin-orbit effects typically appear as quantum cor- 
rections, no sign of spin-orbit coupling can be seen in 
such a calculation. In particular, the magnetic field de- 
pendence of the perturbative density of states is the same 
with and without spin-orbit coupling — the only differ- 
ence being that in a magnetic field, in the absence of 
spin-rotation symmetry, there is no level degeneracy thus 
A gap is 2 2 / 3 times smaller than in the spin-rotation sym- 
metric case. 

To capture the spectral properties on the mesoscopic 
energy scale of order A gap , one needs to go beyond per- 
turbation theory. According to the universality hypoth- 
esis of Vavilov et al*£ the probability distribution of the 
lowest level E\ in properly scaled units is universal and 
identical to the distribution of the smallest eigenvalue 
of random Gaussian matrices from the three symmetry 
classes of RMT. The appropriate scaling is in terms of 
the dimensionless variable x\ — (E\ — E mv ) / A gap and 
the universal distributions are given by^^ 
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where 
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Fi(ar) = ^/F 2 {x)exp J q(x') dx'^j , (5a) 
F 2 (x) = exp (- J (x- x')q(x') 2 dx'^j , (5b) 
(^3) = V^cosh Q J* q(x>) dx'^j . (5c) 



The function q(x) is the solution of the differential equa- 
tion 

q"{x) = -xq{x) +2q 3 (x). 

The boundary condition is q(x) — > Ai(— x) as x — > —00, 
with Ai(x) being the Airy function. The three distribu- 
tions are plotted in Fig. [5] (top panel). The symmetry 
index (3 takes values (3=1 for time-reversal and spin- 
rotation invariant systems, (3 — 4 when time-reversal 
symmetry is present but spin-rotation symmetry is bro- 
ken, and (3 = 2 for systems with broken time-reversal 
symmetry. 

Near the gap, the average density of states in terms of 
the variable x = (E — E gap ) / A gap is given b y 13 i 14 



pi(x) = p 2 (x) + |Ai(-a;) 
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The distribution P 2 and the density of states p 2 
are applicable in an intermediate magnetic field range, 
which exists because the flux needed to close the 
gap is much larger than the flux needed to break 
the time-reversal symmetry. For intermediate fluxes 
> {h/e)N 1 / 3 y / T cl - g 6/h with r erg = V~A/v-p the ergodic 
time and v-p the Fermi velocity), there will still be a gap, 
but its fluctuations are governed by the (3 = 2 symmetry 
class. In this case the presence or absence of spin-orbit 
coupling only affects the parameter A gap (which is re- 
duced by a factor 2 2 / 3 in the absence of spin-orbit cou- 
pling, because the level degeneracy parameter goes from 
s = 2 to s = 1); the gap distribution P 2 and the density 
of states p 2 in rescaled variables do not depend on the 
presence or absence of spin-rotation symmetry. 

If the flux is made much larger ($ 3> 
(h/e)yj Nr CIg S /h) , such that the gap closes, spin- 
orbit coupling starts to play a role again. The reason 
is that an energy level E and its mirror level at — E 
can repel each other, and this repulsion depends on the 
presence or absence of spin-orbit coupling. When there 
is still a gap these levels are widely separated and this 
repulsion is not effectived 

According to Altland and Zirnbauer^ the RMT of an 
Andreev billiard in strong magnetic field is in a new sym- 
metry class called C (D) in the absence (presence) of spin- 
orbit coupling. The average density of states in these 
symmetry classes is 
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The minus (plus) sign should be taken for symmetry class 
C (D). This result expresses the fact that in the absence 
of spin-orbit coupling, a level and its mirror level repel 
each other leading to a vanishing density of states at E = 
0. In the presence of spin-orbit coupling the repulsion 
disappears and levels pile up at the Fermi level, leading 
to a peak in the density of states at E = 0. 



III. SPIN ANDREEV MAP 

To verify these predictions of RMT in a dynamical 
model we combine the general construction of an Andreev 
map 15 with the spin kicked rotato n 16 ' 17 The starting 
point of our discussion is the spin generalized Bogoliubov- 
De Gennes Hamiltonian 18 
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Here H is the single particle Hamiltonian, Ep is the Fermi 
energy, and A is the superconducting pair-potential. The 
operator T stands for time-reversal, and will be specified 
later. 

With the Bogoliubov-De Gennes Hamiltonian as a 
guide we construct the spin Andreev map. First we 
note that if an electron in the normal metal evolves with 
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time-evolution operator F(t), the hole evolves with the 
transformed time-evolution operator TJ 7 '(i)T _1 . Second, 
since we are interested in low energy phenomena, only 
the dynamics on long time scales is important. On time 
scales much larger than T crg , the dynamics can be de- 
scribed as a mapping on a two-dimensional Poincare sur- 
face of section. This amounts to a stroboscopic descrip- 
tion where we are only concerned with the state of the 
electron when it bounces off the boundary. 

The quantum map we use is the computationally ef- 
ficient spin kicked rotator, given in terms of a Floquct 
matrix ; 16 ' 17 

F W = e iso (UUXU^U) w , I, I' = 0, 1, . . . , M - 1. (9) 

The integer M sets the level spacing 5 = 2tt/M. The 
M x M matrices appearing in Eq. @ have quaternion 
matrix elements, and are given by 



Snie -i(M/4sr)VVi*l/M) 
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with 



V{6) = Kcos( 



) (To + K so {(ii sin 20 + <7 3 sin0). (11) 



The quaternions are represented using the Pauli matri- 
ces <Ji with cr the 2x2 unit matrix. The matrix X 
corresponds to the spin-orbit coupled free motion inside 
the dot and n gives scattering off the boundaries of the 
dot. This map is classically chaotic for kicking strength 
K > 7.5. The parameter K so breaks spin-rotation sym- 
metry, 8q breaks time-reversal symmetry and Iq breaks 
other symmetries of the map. The spin-orbit coupling 
time t so (in units of the stroboscopic period tq w T crg ) 
is related to K so through t so = 32n 2 / (K so M) 2 . The pa- 
rameter £o corresponds to the Fermi energy. In the above 
representation of the Floquet matrix, the time-reversal 
operator is given by T = ia%K, where K. is the operator of 
complex conjugation^ Therefore, the hole Floquet ma- 
trix is given by o^F*^ = F, where the overbar denotes 
quaternionic complex conjugation. 

The spin Andreev map is constructed from the electron 
and hole Floquet matrices in the same way as in the 
absence of spin-orbit coupling^ 
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The projection matrix P projects onto the contact with 
the superconductor. Its matrix elements are Pa = 
SmTq Si=i where the set of indices {rii} corresponds 
to the modes coupled to the superconductor. The dwell 
time is therefore Tdwoii = M/N. The corresponding 
Thouless energy is Et — NS/Att — (2rd wc n)~ 1 . As shown 
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Figure 2: Probability distribution of the rescaled excitation 
gap x — (Ei — _B g ap)/Ag ap . Smooth curves (upper panel) 
are the predictions of random-matrix theory, for the three 
symmetry classes j3 — 1,2,4. Histograms are the results of 
the numerical simulation of the spin Andreev map, without 
spin-orbit coupling (zero magnetic field, /3 = 1) and with spin- 
orbit coupling (zero magnetic field, /3 = 4; weak magnetic 
field, /3 = 2). The histograms in the lower panel are plotted 
without any fitting, while in the upper panel -E gap and A gap 
are treated as fit parameters. 



in Ref. ,17], the magnetic field scale at which the gap closes 
is given by 6 C = iiry/N '/(KM 3 / 2 ). From the definitions 
in Eq. ^ the scaling parameters in the spin Andreev 
map become 
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The Floquet matrix has the symmetry 



T = CTF(CT)- 
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corresponding to the CT-antisymmetry of Ti-BdGi the fun- 
damental discrete symmetry of normal-superconducting 
systems £ The eigenphases of the Floquet matrix J 7 , de- 
fined as the solutions of 
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Figure 3: Average density of states in rescaled energy units, 
x = (E — _Z5g a p)/Ag ap , in zero magnetic field (upper panel, 
[3 — 1,4) and in a weak magnetic field (lower panel, = 2). 
The smooth curves are the RMT predictions, the histograms 
are the results of the simulation (using the same values of the 
fit parameters -E gap , A gap as in the upper panel of Fig. [2j. 



play the role of the discrete excitation energies in the 
Andreev billiard. From the symmetry (|14p it follows that 
they come in pairs, ±e, as required. 



IV. NUMERICAL RESULTS AND 
COMPARISON WITH RMT 

In Fig. [2] we plot the excitation gap distribution (his- 
tograms) from our numerical simulation with parame- 
ters K = 41.123, M = 4096, N = 205. The smallest 
e solving Eq. (|15p was calculated for some 6000 differ- 
ent Fermi energies and positions of the contact to the 
superconductor. To generate the three symmetry classes 
we took: 9 /9 c = 0, r dweU /r so = (fi = 1); 9 /9 c = 0, 
Tdweii/T S o = 625 08 = 4); 9 /9 c = 0.4, r dwoU /r so = 625 
(f3 = 2). The values of the parameters c and d for 
9 /9 c = are given by Eq. fl3). For 9 /9 c = 0.4 we calcu- 
late c = 0.427, d = 0.339 from Ref.0. The data is shown 
without any fit parameter in the lower panel and with 
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Figure 4: Average density of states in a magnetic field suffi- 
ciently strong to close the excitation gap. Upper panel: with- 
out spin-orbit coupling; Lower panel: with spin-orbit cou- 
pling. The smooth curves are the RMT predictions from 
Eq. (0 and the histograms are the results of the simulation 
(now without any fit parameters). 



E gap , Ag a p as fit parameters in the u ppe r panel. (A sim- 
ilar fitting procedure was used in Ref.Qjl) The fitted val- 
ues of -Egap and A gap do not vary much (by about 5% and 
10%, respectively) from their nominal values [given by 
Eq. (|13p ], but the agreement with the RMT predictions 
improves considerably if we allow for this variation^ 

The characteristics of the gap distribution are clearly 
obtained in our simulation. In zero magnetic field, the 
gap distribution becomes narrower as the strength of 
spin-orbit coupling is increased (0 = 1 — > 4). The RMT 
prediction for the standard deviation of the scaled dis- 
tributions Pp(x) is o\ — 1.27 and o± = 0.64. The cor- 
responding values obtained in our numerical simulation 
without any fitting are <j\ = 1.34 and (74 = 0.72. If a weak 
magnetic field is present ((3 — 2), the width of the distri- 
bution is predicted to be intermediate between the cases 
with j3 = 1,4. As seen, the dynamical model follows the 
prediction. The theoretical value of the standard devia- 
tion is o~2 = 0.90, the numerical result (without fitting) 
is ct 2 = 0.99. 

Using the same fit parameters as in Fig. [2] we plot the 
average density of states close to the gap in Fig. [3J As 
seen, the numerical data follow closely the analytical pre- 
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dictions, the deviations becoming significant only outside 
the universal regime \E — E ga , p \ -C E T , i.e. |x| <C N 2 / 3 . 
In the absence of magnetic field the spin-orbit coupling 
induced oscillations are clearly obtained. 

The density of states in a strong magnetic field is given 
in Fig. 21 The upper panel shows the data without spin- 
orbit coupling, 5 and the lower panel shows what happens 
if spin-rotation symmetry is broken. The numerical data 
follows closely the analytical prediction (J7J) of RMT. In 
particular the enhanced density of states in the presence 
of spin-orbit coupling is clearly seen. The first oscillations 
in the density of states are also captured in the dynam- 
ical model. The frequency doubling due to the reduced 
degeneracy (s = 2 — > s — 1) is apparent. 

V. CONCLUSION 

In conclusion, we have introduced a quantum map 
for the dynamics of a chaotic quantum dot with spin- 
orbit coupling connected to a superconductor. We have 
demonstrated three effects of spin-orbit coupling on the 
excitation spectrum of this Andreev billiard: The nar- 



rowing of the distribution of the excitation gap and the 
appearance of oscillations in the density of states in the 
absence of a magnetic field; and the peak in the density 
of states at the Fermi level in strong magnetic field. Our 
numerical simulations confirm the predictions of random- 
matrix theory. The third effect is particularly interest- 
ing from an experimental point of view. In view of the 
possibility to tune the strength of spin-orbit coupling in 
quantum dots ; 21 ' 22 one can imagine tuning the density 
of states at the Fermi level from zero to a value of twice 
the normal density of states. 
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